The Status of Black-Hole Binary Merger Simulations with Numerical Relativity 
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The advent of long-term stability in numerical relativity has yielded a windfall of answers to 
long-standing questions regarding the dynamics of space-time, matter, and electromagnetic fields in 
the strong-field regime of black-hole binary mergers. In this review, we will briefly summarize the 
methodology currently applied to these problems, emphasizing the most recent advancements. We 
will discuss recent results of astrophysical relevance, and present some novel interpretation. Though 
we primarily present a review, we also present a simple analytical model for the time-dependent 
Poynting flux from two orbiting black holes immersed in a magnetic field, which compares favorably 
with recent numerical results. Finally, we will discuss recent advancements in our theoretical un- 
derstanding of merger dynamics and gravitational waveforms that have resulted from interpreting 
the ever-growing body of numerical relativity results. 
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I. INTRODUCTION 

In the last five years, the field of gravitational wave 
(GW) astronomy has witnessed a remarkable conver- 
gence of experimental and theoretical achievement. The 
current generation of GW interferometers, particularly 
the Laser Interferometer Gravitational Wave Observa- 
tory (LIGO) and VIRGO detectors, has positioned this 
community on the cusp of discovery, by achieving their 
design sensitivities and making the first direct detection 
of GWs a real possibility. Meanwhile, the GW signature 
from the merger of a black- hole binary (BHB), expected 
to be the most common source for these and future ob- 
servatories, went from being a complete unknown to a 
reasonably well-understood, surprisingly smooth transi- 
tion between the earlier inspiral and the final ringdown, 
thanks to tremendous advancements in the field of nu- 
merical relativity (NR). In this article, we will limit our 
scope to the advancements that have occurred in just the 
past year, with only a minimal amount of background 
information, in order to provide a snapshot of the most 
current research in the field. For a review of progress in 
the preceding year, we refer the reader to [l|, and for a 
broader review of NR achievements, please see 0, Q. We 
consider it an indicator of the healthy progress of the field 
that, in addition to the general review articles typical of 
any field, in recent years the field of NR has progressed 
rapidly enough to warrant an annual review of the past 
year's progress. 

In Sec. ini we give a brief review of the equations be- 
ing solved by numerical relativists, and the methods that 
are most typically applied to solve them. In Sec. Illli we 
discuss the most recent achievements in NR methodol- 
ogy, and the simulations of numerically challenging sys- 
tems of astrophysical interest that have been facilitated. 



In Sec. IIVI we briefly mention novel theoretical devel- 
opments resulting from the interpretation of numerical 
results. We summarize and conclude in Sec. |Vl 



II. STANDARD METHODOLOGY 

The community of numerical relativists attempts to 
solve Einstein's equations in the strong-field regime 
numerically (hence the name). During the early, 
approximately-adiabatic inspiral, we can treat the sys- 
tem perturbatively through expansions in v/c or other 
small parameters related to it. After the merger, we can 
again treat the system perturbatively, as the quasinor- 
mal ringing of small deviations from a Kerr background. 
However, during the merger, unless the masses are dis- 
parate and q = mi/m2 provides a tractably small ex- 
pansion parameter, we can no longer treat the system 
perturbatively, and are forced to solve the full equations, 
given by 



(1) 
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where R^i,, R, g^^, and T^^ are the Ricci tensor, its trace, 
and the metric and stress-energy tensors, respectively. 

In the case of BHBs in vacuum, T^j^ vanishes. The 
resulting equation is generally evolved by splitting up 
space-time into a sequence of three-dimensional hyper- 
surfaces. The choice of how to make this split is not 
unique, as different evolution equations with different 
variables can be used, as well as different coordinates. 
The coordinate freedom is fixed through the choice of 
dividing space-time into hypersurfaces. The lapse, a, 
determines the coordinate time interval between one hy- 
persurface and the next. The shift, f3i, determines the 
motion of the three spatial coordinates within each hy- 
persurface. The splitting of space-time results in a set of 
evolution equations, plus constraints that need to be sat- 
isfied. To the extent that the constraints vanish, all evo- 
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lution equations must be analytically equivalent. How- 
ever, different choices will have different numerical prop- 
erties. Indeed, the four decades that passed between the 
first attempts to perform NR simulations and the recent 
breakthroughs were spent, in large part, searching for 
formulations that would be numerically stable. 

At present, there are primarily three analytic forms 
for Einstein's equations that are actively being applied 
to the BHB problem: the BSSN formulation [iBand 
two different Generalized Harmonic formulations [S-Q- 
Most typically, the equations are spatially discretized 
with finite differences (with the exception of SpEC [ioI |. 
which uses pseudo-spectral methods). The state-of-the- 
art codes employ eighth-order spatial differencing, but 
most codes have a mixture of orders for different parts 
of the code. All implementations of which we are aware 
discretize in time using Runge-Kutta methods, usually 
to fourth-order accuracy in the time step. Given the 
disparate scales involved in the problem (r > lOOAf for 
the wave zone, A < M/lOO for the near- field), a range of 
resolutions within each simulation is required, and this 
is achieved in a number of different ways by different 
groups. There are man y c odes being actively used and 
developed 0, Hi, [H, [13, [Hi, [i3> and 

comparisons have been published that verify the consis- 
tency of results among different groups [l8lfl9|. 

In the presence of matter and/or electromagnetic (EM) 
fields, the equations governing the evolution of space- 
time must be coupled to an appropriate stress-energy 
tensor. It is noteworthy that matter is not needed for 
T^jy to be nonvanishing. In the case of the Einstein- 
Maxwell system (see e. g. [13), the sources for EM fields 
are not included in the stress-energy tensor. The in- 
clusion of these sources, the fields themselves, and the 
equations of hydrodynamics, all coupled to Einstein's 
equations, constitute the equations of general relativistic 
magneto- hydrodynamics (GRMHD), while the coupling 
of the hydrodynamics and Einstein's equations, without 
EM fields, make up the general relativistic hydrodynam- 
ics (GRHD) equations. The recent advent of GR(M)HD 
codes has yielded some of the most interesting science in 
NR, and perhaps in all of theoretical astrophysics, in the 
past year, and will be discussed in Sec. IIIIDI 



III. NOVEL METHODOLOGY AND 
ASTROPHYSICS 

In this section, we will discuss novel developments in 
NR methodology, and the novel astrophysical results that 
those developments have facilitated. 



perturbative methods, and they are, at present, com- 
pletely impossible to simulate non-perturbatively. How- 
ever, for more intermediate mass ratios, (7> 1/100, it is 
now possible to perform short simulations which include 
1-2 orbits of the late inspiral, the merger, and the ring- 
down. This is an impressive achievement, because the 
resolution demands are set by the smaller black hole, so 
that a black hole 100 times smaller than its companion 
requires 100 times finer resolution. The trailblazers in 
the endeavor to simulate these systems have been the 
RIT group, using their LazEv code [2l|. Their first ef- 
fort involved combining a full numerical evolution of the 
punctures with a perturbative treatment of the radiation 
[22l |. Later, they validated this approach by comparing 
with fully numerical simulations for a g = 1/10 system 
[13. Most recently, they have achieved the aforemen- 
tioned 1/100 mass ratio through a combination of im- 
proved gauge conditions, an improved mesh refinement 
scheme, and an improved allotment of cpu cycles [24l|. 



B. Toward Maximal Kerr Initial Data 

The Caltech-Cornell-CITA collaboration has overcome 
the obstacles of generically stable black-hole binary evo- 
lutions, and have demonstrated the ability to stably 
evolve systems with moderate mass ratios q<l/2 (25| 
and extremely high spins a = 5/M^<0.95 [13, in- 
cluding merger and ringdown, using their SpEC code 
[13 . To facilitate the simulation of stable near-extremal 
spins, they had to employ novel methodology in order 
to exceed the bound of a = S/M'^< 0.93 that IS an 
intrinsic characteristic of Bowen-York initial data 27]. 
More specifically, while the spin in Bowen-York initial 
data can be set arbitrarily high, it will quickly relax to 
the aforementioned bound, so that larger spins cannot 
be evolved with that approach. To this end, the au- 
thors of [23 applied a novel method [13 wherein the 
extended-confornial-thin-sandwich equations [13 were 
combined with a conformally-curved metric, constructed 
through the weighted superposition of the metrics for two 
boosted, spinning black holes in Kerr-Schild coordinates. 
Using this approach facilitated the simulation of a pair of 
equal-mass black holes with spins of a = 0.95 anti-aligned 
with the orbital angular momentum (thereby minimizing 
the total system angular momentum) [26|. We note that 
the Illinois group also developed a novel approach for 
near-extremal initial data [33 . but it has not yet been 
demonstrated in a binary evolution. 



C. Nonlinear memory 



A. "Extreme" mass ratios 

Whereas the typical definition of extreme mass ratio 
systems is q ~ 10~^, such systems are of little interest to 
numerical relativists, as they can be treated well with 



An interesting consequence of the nonlinear nature of 
GW propagation is the emission of GWs generated by 
GWs. This effect has long been established analytically, 
and is often referred to as the "Christodoulou memory" 
[3l| . However, this effect is exceedingly small, and is 
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therefore difScult to simulate, as systematic error in the 
simulations must be smaller than the effect. The GW 
memory, though small, could conceivably be detectable 
by the space-based Laser Interferometer Space Antenna 
(LISA) [32|- The authors of ^] were able to accurately 
calculate the memory effect from the final merger and 
ringdown, by using novel numerical methodology imple- 
mented in their Llama code, which includes multipatch 
capabilities fisl . [sj] and Cauchy characteristic extraction 
(CCE) [H, H^. The multipatch code facilitates multi- 
ple coordinate patches in different regions, allowing tra- 
ditional Cartesian grids around the near zone, with a 
smooth transition to a six-patch polar grid in the wave 
zone. Since the waves do not gain angular structure in 
the wavezone, this allows one to increase the grid size 
by only adding radial elements, thereby greatly decreas- 
ing the cost of large grids, and permitting much more 
accuracy in the region where the waves propagate. The 
implementation of CCE represents the first fully gauge 
invariant GW calculation, wherein the GWs are formally 
measured at infinity, rather than extrapolating to infin- 
ity from a sequence of finite extraction radii (though we 
note that progress has also been made in the last year on 
compactification, including infinity explicitly in the com- 
putational domain [131). Using Llama, the authors of 
[33| were able to accurately resolve the nonlinear memory 
for the first time. We note that, in [H, 111], the /i2o com- 
ponent of the GW strain (which contains the memory) 
appears to be a significant fraction of the fundamental 
quadrupole modes, so that some readers might be con- 
fused why the signal-to-noise ratio (SNR) of the memory 
contribution is roughly two orders of magnitude below 
the SNR of the total signal (see Fig. 2 of [3l|). In this 
context, one needs to bear in mind that constant offsets 
in strain are gauge-dependent and undetectable. There- 
fore, only the derivative of the memory contribution to 
the GW strain contributes to the SNR [H. 



D. GR(M)HD 

Perhaps the most exciting developments have been the 
coming-of-age of GR(M)HD codes. There are currently 
four codes capable of GRMHD with dynamical space- 
time: WhiskyMHD [11], SACRA [13, the LSU-LIU- 
BYU-PI collaboration code [sl], and the Illinois group 
code [i^. A number of other groups have also devel- 
oped, or are currently developing, working GRHD codes. 
As many results of astrophysical interest do not require 
including electromagnetism (EM), these codes are also 
an exciting development. In addition to implementing 
several stable evolution schemes that satisfy the con- 
straints of Einstein's equations along with other consis- 
tency requirements (such as preserving the divergence- 
free character of the magnetic field for GRMHD), most 
of the aforementioned codes have either incorporated, 
or are currently developing, more relevant physics for 
problems involving white dwarfs, neutron stars, and any 



gaseous environment, such as photon and neutrino trans- 
port, and stellar equations-of-state (EOS), ranging from 
simple polytropes to microphysical EOSs. While most 
of the work in this field has so far centered on mergers 
of binaries that include a neutron star or white dwarf, 
and is therefore beyond the scope of this review, there 
have been a few results from GR(M)IID codes applied to 
BHBs in gaseous environments. 

The Illinois group estimated the accretion luminosity 
from a BHB immersed is a gaseous environment using 
GRHD, and compared the result with the luminosity of 
a single black hole of the same total mass, immersed in 
the same environment [4l| . By estimating the contri- 
butions from bremsstrahlung and synchrotron emission, 
they predicted an enhancement of luminosity by three 
orders of magnitude at BHB merger, compared to the 
single-black-hole case. If true, this result would have a 
significant impact on the search for EM counterparts of 
GW signals from BHBs. 

In addition to accretion luminosity from a surrounding 
disk, significant luminosity may accompany BHB merg- 
ers in the form of jets, particularly in the presence of 
an EM field. Much progress has been made in studying 
this possibility. In [201], which built upon earlier studies 
[39t, ,42j , the Einstein-Maxwell equations were evolved for 
the case of a BHB immersed in an EM field, with differ- 
ent BHB spin configurations. By calculating the Poynt- 
ing flux, they showed that the direct EM emission was 13 
orders of magnitude smaller than the GW emission. Fur- 
thermore, as the EM frequency evolution tracked that for 
GWs, the direct EM signature would occur in the mHz 
range for supermassive BHBs, and would therefore be at 
far too low a frequency to be detectable. 

However, by combining the evolution of the Einstein- 
Maxwell equations with a tenuous plasma, which is 
evolved using the force-free approximation [isl |43 ] 
wherein the inertia of the plasma is neglected, the authors 
of [i^ HI] were able to calculate the synchrotron radia- 
tion resulting from the acceleration of the plasma due to 
the EM Poynting flux previously discussed, upconverting 
that energy to GHz, and thereby making it detectable by 
X-ray observatories. The radiation transitions from an 
TO = 2 multipolar structure to an to = structure, which 
has implications for determining whether the EM emis- 
sion is in fact of the same origin as a given GW observa- 
tion, since there will likely be many variable EM sources 
within the error ellipse of any GW observation. How- 
ever, the most interesting aspect of [45|, from a theoreti- 
cal standpoint, is the fact that the EM emission is gener- 
ated by a BHB consisting of two nonspinning black holes. 
Since the emission from a single black hole is most typi- 
cally associated with the Blandford-Znajek (BZ) mecha- 
nism (43| , which requires that the black hole be spinning 
in order to operate, this new result may be surprising at 
first. The flux in the BZ process originates from the fact 
that a spinning black hole will twist magnetic field lines, 
which are assumed to be anchored to a distant accretion 
disk. This twisting gives rise to the EM Poynting flux. 
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with the torquing of the magnetic field lines causing a 
decrease in the black hole's spin. 

In the case of a BHB, it has now been shown that 
there is a significant contribution to EM emission that is 
a generic feature, depending only on the presence of an 
orbiting binary and an EM field. Indeed, given the obser- 
vation of protostellar jets and their apparent similarities 
to their extragalactic counterparts 47] , it is possible that 
the conditions for jet formation may be far more generic 
than was previously thought, with the presence of plasma 
and an appropriate magnetic field strength and geometry 
being the only clear prerequisites. 

In the context of BHBs, the existence of jets, indepen- 
dent of spin, may not be as surprising when thought of 
within the framework of the membrane paradigm JS], 
wherein the black hole event horizon is reinterpreted as 
a two-dimensional viscous membrane, and the standard 
Maxwell equations in three dimensions govern the dy- 
namics of EM fields evolving along a sequence of hyper- 
surfaces. It is well understood that, from this viewpoint, 
the BZ effect can be calculated using simple circuit equa- 
tions. In this picture, a magnetic field threading the black 
hole and frozen in to a distant accretion disk serves the 
role of electrical wires, with charges spiraling along them. 
The disk and the horizon membrane serve as resistors, 
and the potential results from the twisting of the mag- 
netic field lines by the spin of the black hole. 

The authors of (1^, Hg] suggested that, in the case of 
a BHB, the Faraday induction results from the orbital 
motion of the two black holes through the fixed magnetic 
field, and could also be understood through application of 
the membrane paradigm. We carry out the calculation 
here, which does indeed bear out the scaling behavior 
suggested in [4^, . Assuming that the plasma is low- 
density and the force-free approximation applies, that 
the orbit is in the {p, 0} plane with orbital velocity v 
(though any relative velocity between the black hole and 
the magnetic field will have the same effect), and that 
the background magnetic field is given by -B = Bz (in 
cylindrical coordinates {p, 0, z}), the induced potential 
from one orbiting hole is then given by 



y = ^ adl- E = - <j> adl- {v x B) = -2aruvB 



(2) 



where a = ^1 - 2M / r is the lapse between hypersur- 
faces of constant Schwarzschild time, M the combined 
rest mass of both black holes, and rn = M is the horizon 
radius of either hole. The contour integration for employ- 
ing Faraday's law in Eq. [2] is somewhat subtle. We place 
one side along one hemisphere of the stretched horizon, 
and assume the opposing side resides far away, where the 
B-field is weak, so that its contribution to the integral 
vanishes. We note that this region is assumed to be far 
away from the computational domain in [4^, |4^ , where 
the B-field is uniform. The other sides of the contour 
are perpendicular to v x B where it remains uniform, so 
that those contributions also vanish, and only the inte- 
gral along the horizon remains. 



To leading order, the orbital velocity is given by v 
^ M/r, so that Eq. [2] can be re-expressed as 



V = -2BM 




(3) 



Using the remarkable result that the membrane has an 
effective resistance, i?H, of 377 (or 47r in geometrized 
units) '49'], and bearing in mind that each hole emits in 
both hemispheres, we can solve for the total Poynting 
flux from both orbiting holes: 



4 o ,Af 



I 



(4) 



We note that it is often assumed that the astrophysical 
load also contributes a resistance, Rl ~ Rh [H,!!!], but 
as this is rather ad hoc and depends on poorly understood 
astrophysics, and there is no apparent way that a distant 
astrophysical load could manifest itself in the simulations 
of [4^, |4y] , we neglect this contribution. Eq. 2] predicts a 
peak luminosity at r = 4M, which we can express in cgs 
units with appropriate scaling relationships as 



I 

2^^ 



BZ 



2x10*-^ erg/s 



B 
104 G 



M 



108 M, 



, (5) 



where a is the dimensionless spin parameter of the Kerr 
black hole for the single black hole BZ process (0 < a < 1), 
and Lbz is understood to be the luminosity from the 
BZ process for a single black hole of spin a and mass 
M (to approximate the mass of the merged remnant). 
We emphasize that this estimate is only valid prior to 
merger, and indeed predicts a vanishing luminosity for 
r = 2Af, so that the merger and its accompanying burst 
of luminosity seen in [45l.T46j are not included. However, 
if we combine Eq. with the leading order relationship 
for r{t), 

we see excellent agreement between our simple model and 
the data from Fig. 4 in [i^ (see Fig. [T]). 

We also note that this model differs significantly from 
that suggested in [s^]- Using Eq. 7 in 50] to find the 
prediction for Lmax in that model, we find that the pre- 
diction is approximately four orders of magnitude larger 
than the numerical result found in (45i] . While neither 
model should be expected to perform well in predicting 
the true peak luminosity, which occurs in a highly dy- 
namical and strong-field regime, the suggestion in [50] 
that L cx M^B^ yields vastly different predictions from 
Eq. m and disagrees with the scaling suggested in [4^ . 



E. State of the art 

It is useful to summarize the boundaries of parameter 
space that have currently been explored. While some- 
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FIG. 1: Comparison of the numerical result from Fig. 4 of 
[4^ (solid line) for the Poynting flux from two orbiting non- 
spinning black holes with the simple model given by Eqs. |4] 
and [6] (dashed line). 



Most waveform cycles 


64 


[51] 


Smallest mass ratio 


1/100 


[24] 


Largest kick (astrophysical) 


3254 km/s 


[52] 


Largest kick (hyperbolic) 


9589 km/s 


[53] 


Largest initial spin (So/M^dm) 


0.97 


[28] 


Largest final spin {Si/Mj^j2>M) 


0.96 


[54] 


Lowest eccentricity 


0(10"^) 


[55] 


Most energetic (-Erad/MAOM) 


0.35 


[54] 



TABLE I: Summary of the most extreme BHB systems that 
have been evolved to-date in numerical relativity simulations, 
including the extreme quantity, its value, and the reference 
for that simulation. 



what arbitrary, a summary of the extreme cases of dif- 
ferent quantities that have been simulated is useful and 
serves as a snapshot of the current status of the field. We 
include a set of physical quantities from NR simulations, 
the superlative result for that quantity, and the reference 
wherein that result was found, in Table ID 

The fact that a number of the extremes have not oc- 
curred within the past year is perhaps indicative that 
the tide of NR breakthroughs that began in 2005 has 
begun to ebb. However, as the pace of 2005-2009 could 
not be maintained indefinitely, it is our hope that the 
slowdown is simply the start of a period of more mod- 
est, but sustainable, growth in the field. We emphasize 
that the list of superlatives is restricted to results that 
have been achieved in full BHB simulations, as opposed 
to results from single-black-hole simulations, or extrapo- 
lations of a set of BHB simulations to some extreme. We 
also emphasize that we are summarizing the extremes in 
physical systems that have been simulated, and not in- 
cluding extremes in the numerics (such as smallest phase 
error, etc.), as any set of such results is arbitrary, and 
difficult to summarize in an unambiguous way. 



IV. THEORETICAL ADVANCES 

Apart from achievements in simulating new systems, a 
burgeoning subfield of NR is the ongoing effort to draw 
general theoretical insights from the ever-growing body of 
available NR simulations. Independent research projects 
at Caltech [H, 113] and Chicago [llj have endeavored to 
gain analytical insight into the nonlinear dynamics of mo- 
mentum transfer in merging BHBs, particularly with re- 
spect to the bobbing and ultimate recoil. Substantial 
work has been done to condense the body of simulations 
into a simple fitting formula, based on either the final 
state of the binary just prior to merger j59| or on the 
initial conditions at wider separations [60[. More gen- 
eral efforts to understand the complete final state of the 
merged remnant, including the final mass and spin, have 
been successful in predicting subsequent numerical re- 
sults [HI,!!!. 

Progress has also been made in moving beyond purely 
phenomenological fits of merger waveforms. Previously, 
progress was made by fitting to natural extensions of 
the post-Newtonian |63| or effective-one-body formalisms 
[slilllj. More recently, attempts have been made to ap- 
ply our improved theoretical understanding of the dom- 
inant dynamics during merger in order to formulate 
more physically-motivated merger waveforms. These at- 
tempts are an effort to move beyond simply bridging the 
gap between the well-understood inspiral and the well- 
understood ringdown, in order to encapsulate the key 
features of generic mergers, in the hopes that the results 
will prove to be predictive for as-yet-unsimulated (or cur- 
rently unsimulatable) systems. In [S^,!!^], the close-limit 
approximation was applied in an effort to predict the re- 
coil of comparable-mass BHBs. This research program 
has a long history, but the novelty in these investiga- 
tions is the quality and analyticity of the post-Newtonian 
initial data, and the existence of full NR results, which 
largely validate the effectiveness of this procedure. Simi- 
larly, in [i^, a hybrid of post- Newtonian and black hole 
perturbation theory techniques was used to study the 
radiated energy and momentum and the gravitational 
waveforms for head-on collisions. The results compared 
favorably with full NR simulations. 

In [i^ [t^I, the "impficit rotating source" model for 
merger waveforms was developed, and its results were 
rigorously compared with full numerical simulations for 
nonspinning systems with mass ratios g<l/6. This 
model was built on the observation that merger wave- 
forms from numerical simulations behave as though they 
were generated by a rigid rotator. Specifically, the fre- 
quency evolution of all £ = m modes, when weighted by 
m, are equal across all modes and mass ratios, consistent 
with multipoles of a single object rotating with a well- 
defined frequency. The angular momentum at merger is 
also proportional to the frequency, as one would expect 
for a rotator with a well-defined moment of inertia. These 
and other indicators led the authors of to formulate a 
physically-motivated merger waveform, which naturally 
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ties to the final ringdown phase of the waveform, whereas 
most other methods tie to the late inspiral waveform, at a 
point where it may have already accumulated a large de- 
gree of error. Combined with the aforementioned models 
for the final mass and spin of the merged remnant across 
parameter space, this model may provide a novel method 
for predicting the merger waveform for generic systems, 
with a minimum of fine-tuning. 

The commonality of all the approaches to gaining an- 
alytical insight into merger dynamics is the observation 
that the merger, though happening in the strong-field, 
nonlinear regime, seems to be well-described by tools de- 
veloped primarily for linear pcrturbative analyses in the 
weak-field regime. This observation is one of the most 
fundamental results that has been made possible by nu- 
merical relativity. 



significant astrophysical and theoretical interest. In the 
past year, progress has been made in simulating longer 
waveforms overall, as well as systems with more extreme 
mass ratios and larger spins. Advances in methodology 
have significantly enhanced the accuracy and efficiency 
achievable by state-of-the-art codes. Progress has also 
been made in using the body of available numerical sim- 
ulations to inform a greater theoretical understanding 
of the strong-field behavior of black-hole binaries. The 
inclusion of matter and electromagnetic fields in black- 
hole binary simulations has led to unexpected discoveries, 
which will have significant astrophysical implications in 
the years to come. The field of numerical relativity has 
now entered a period of more modest, but hopefully sus- 
tainable growth, with a substantial amount of discovery 
space that remains to be explored. 



V. CONCLUSIONS 

We have presented a brief overview of the current 
state of black-hole binary simulations in numerical rel- 
ativity. While progress has slowed somewhat since the 
"gold rush" period following the breakthroughs in 2005- 
2006, there have still been a number of novel results of 
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